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We analyze theoretically the experiment reported in [F. Gerbier et al, )coml-mat/0503452l , where the in- 
terference pattern produced by an expanding atomic cloud in the Mott insulator regime was observed. This 
interference pattern, indicative of short-range coherence in the system, could be traced back to the presence 
of a small amount of particle/hole pairs in the insulating phase for finite lattice depths. In this paper, we an- 
alyze the influence of these pairs on the interference pattern using a random phase approximation, and derive 
the corresponding visibility. We also account for the inhomogeneity inherent to atom traps in a local density 
approximation. The calculations reproduce the experimental observations, except for very large lattice depths. 
The deviation from the measurement in this range is attributed to the increasing importance of non-adiabatic 
effects. 
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The superfluid to Mott insulator (MI) transition undergone 
by an ultracold Bose gas in an optical lattice has attracted 
much attention in the recent years as a prototype for strongly 
correlated quantum phases Ql, |2. yl 01 ■ A key observable in 
these systems is the interference pattern observed after releas- 
ing the gas from the lattice and letting it expand for a certain 
time of flight. Monitoring the evolution of this interference 
pattern not only reveals the superfluid-to-MI transition Q,0], 
but also allows for example the detection of number-squeezed 
states in the lattice ||5K6|i|, or the observation of collapse and 
revivals of coherence due to atomic interactions 0]. Because 
of its experimental importance, a quantitative understanding 
of this interference signal is crucial to characterize quantum 
phases of bosons in optical lattices. 

Although no interference pattern is expected for a uniform 
array of Fock states (what we call a "perfect" Mott Insula- 
tor) l29il . a finite visibility is nevertheless observed in exper- 
iments above the insulator transition f^HSl, in agreement 
with numerical calculations [,9, 10|. We have studied this phe- 
nomenon experimentally, and shown that despite its insulating 
nature that forbids long-range coherence, a MI still exhibits 
short-range coherence at the scale of a few lattice sites |8]. 
This can be attributed to the structure of the ground state for 
finite lattice depths, which consists of a small admixture of 
particle/hole pairs on top of a perfect MI. A qualitative model 
based on a lowest-order calculation of the ground state wave- 
function was also presented in our previous work, which re- 
produced the main trend and order of magnitude of the ob- 
served visibility. 

In the present paper, we would like to present a more pre- 
cise calculation that includes higher order corrections (see 
also fllll '). We describe a MI state at zero temperature using 
the Random Phase Approximation (RPA), already introduced 
in Refs. (H [H Q QJ d ■ Instead of the path integral ap- 
proach used by these authors, we obtain here the RPA Green's 
function using a diff'erent method inspired by Hubbard's orig- 
inal treatment of the fermionic model Taking the 
experimental geometry and the inhomogeneous particle distri- 
bution into account, we find good agreement with our exper- 
imental data, which provides further support for the physical 
picture presented above. 



The paper is organized as follows. In section |I| we recall 
the description of ultracold atoms in an optical lattice by the 
Bose-Hubbard model, and discuss the inhomogeneous shell 
structure that develops in an external confining potential. Sec- 
tion lUpresents the calculation of the interference pattern ob- 
served after free expansion of the atom cloud and its link with 
the quasi-momentum distribution. The main results are pre- 
sented in sections Hill and [Tvl where we respectively present 
the calculation of the interference pattern in the uniform case 
using the RPA, and extend it to the inhomogeneous case to 
compare to the experimental data of 0. Details of the calcu- 
lation are described in the appendix. 

I. BOSE-HUBBARD HAMILTONIAN 

In this section, we briefly recall the theoretical description 
of an ultracold atomic gas trapped in an optical lattice. The 
optical lattice potential, which results from the superposition 
of three orthogonal and independent pairs of counterpropagat- 
ing laser beams, can be written as 

VoL(r) = Vo (sin2(^Lx) + ^in\ki^y) + sm\k^z)) , (1) 

Here Vo is the lattice depth, k]^ - Itt/A]^ is the laser wavevec- 
tor, /11 is the laser wavelength and m is the atomic mass. 
As usual, we measure Vo in units of the single-photon recoil 
energy Er - h^/2mA^. The lattice potential has a simple 
cubic periodicity in three dimensions, with a lattice spacing 
d = A]^/2 a 425 nm in our case. As shown in Ql, the behavior 
of the atomic system in such a potential can be described by 
the Bose-Hubbard model, defined by the hamiltonian 

9i = -tJ]a'laj + Y,jni(hi-l). (2) 

Here the operator aj creates an atom at site ;, «, = a^ffi, is the 
on-site number operator, and the notation (/, j) restricts the 
sum to nearest neighbors only. The relative strength between 
the tunneling matrix element t and the on-site interaction en- 
ergy U is controlled by the depth Vo of the periodic potential 
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which confines the atoms 13011 . The phase diagram of this 
hamiltonian is well known: The system is in a Ml state within 
characteristic lobes in a f /[/ versus chemical potential fi phase 
diagram, and in a superfluid state outside of these lobes fl9ll . 

In the experiments, an additional potential VextCr) is super- 
imposed to the lattice potential, leading to a spatially varying 
chemical potential across the cloud. This favors the formation 
of a "wedding cake" structure of alternating MI and superfluid 
shells, which reflects the phase diagram of the Bose-Hubbard 
model Jlll^ figlEoll . The external potential is due to a combi- 
nation of a magnetic potential in which the condensate is ini- 
tially formed and of an optical potential due to the Gaussian 
shape of the lattice beams. To a good approximation, it can be 
considered as a harmonic potential with trapping frequency 
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where Wn, is the frequency of the magnetic trap, assumed 
isotropic, and where w is the waist {l/e- radius of the intensity 
profile) of the lattice beams, assumed identical for all axes. 
For large lattice depths, the confinement is mainly due to the 
optical part. 

In the current experiments, this external potential varies 
slowly across the lattice. In this limit, the shell structure 
can be calculated in a local density approximation, which as- 
sumes a known relation «h[yu] between the density n and the 
chemical potential fi for the homogeneous system. Then, the 
coarse-grained density llsill for the inhomogeneous system is 
calculated as «(r) = «h[/U - V'ext(r)]. The chemical potential 
is fixed by the relation = J d'-^^r «(r). For a fixed lattice 
depth and atom number, we calculate numerically the relation 
«h[/i] using mean-field th eory at zero temperature, i.e. in the 
mean-field ground state lll2l 12111 . We then repeat the steps 
outlined above, varying the chemical potential until the target 
atom number is obtained within 0.1 %. For all calculations, 
the values a>m = 27r x 15 Hz and w = 130 fim are used, which 
match our experimental parameters. We show in Fig. an ex- 
ample of such a calculation for a lattice depth of Vq = 18 £r 
and = 2.2 x 10^ atoms. 

The presence of the external potential significantly affects 
the atom distribution in the lattice, which is determined by 
the competition between interaction and potential energy. On 
the one hand, expanding the cloud minimizes the density and 
the interaction energy, and on the other, contracting it mini- 
mizes the potential energy, as in conventional harmonic traps 
1I22I I23I1 . The latter is favored at low atom numbers, where 
only a no = 1 shell forms. When more atoms are added, the 
radius of the unity filled MI region increases until a critical 
atom number, which we estimate to be 7 x lO'* atoms for our 
parameters, is reached. Above this critical atom number, a 
higher density core appears near the trap center A MI with 2 
atoms per site is then obtained near the trap center if the lattice 
depth is above the critical value Vo ~ 14.7 Er. This value has 
been calculated using the boundary derived in 
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where z = 6 is the number of nearest neighbors in 3 dimen- 
sions. In the specific example shown in Fig.^ we have chosen 
Vq - 18 and N -2.2x 10^ atoms, so that both «o - 1 and 
«o = 2 MI are present. Similarly, we calculate for our experi- 
mental parameters that an «o = 3 shell is also present for atom 
numbers larger than 2.7 x 10^, and lattice depths larger than 
16 £r. 
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FIG. 1: Calculated density profile for a lattice depth Vq = IS Eg, and 
A' = 2.2 X 10^ atoms. Here d is the lattice spacing. The upper solid 
line indicates the numerical calculation of the total density, and the 
dotted line is the i7 — » 00 extrapolation (see section lTvl . 



II. INTERFERENCE PATTERN 

We now turn to the description of the interference pattern 
observed after release of the atom cloud from the optical lat- 
tice and a period of free expansion. From an absorption image 
of such a pattern, the phase coherence of the atomic sample 
can be directly probed. The density distribution of the expand- 
ing cloud after a time of flight ff can be calculated as (3013 
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which mirrors the momentum distribution of the original 
cloud. Momentum space and real space in the image plane 
are related by the scaling factor tiujm - independent of the 
lattice parameters. The envelope \w'^ is the Fourier transform 
of the Wannier function in the lowest Bloch band. When each 
potential well is approximated as an harmonic potential, the 
Wannier function is the corresponding gaussian ground state 
wavefunction. The envelope function in Eq. (|5|i then reads 
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where w,f = fitf/mwo, and where wq is the size of the on-site 
Wannier function. Finally, we have defined the quantity 



,/k-(r,-r,) 
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(7) 
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When k is restricted to the first Brillouin zone, S(k) is nothing 
else than the quasi-momentum distribution. Information about 
the many-body system is contained in this quantity, which is 
periodic with the periodicity of the reciprocal lattice In/d. 
Thus, to predict the interference pattern and compare to the 
experiments, our goal is to calculate S(k) for a given lattice 
depth and density. 



III. QUASI-MOMENTUM DISTRIBUTION IN THE 
HOMOGENEOUS MOTT INSULATOR 



even "stretch" over a few lattice sites. This mobility acquired 
by particle/hole pairs is reflected in the modified dispersion re- 
lation (|9|l, which explicitly includes the band structure. Note 
finally that higher order excitations, corresponding to occupa- 
tion numbers no ± 2, no ± 3, ... are neglected. At zero tempera- 
ture, such excitations become important only very close to the 
superfluid transition where the MI is destroyed. 

The quasi-momentum distribution can be di- 
rectly deduced using the general relation S(k) - 
-;limi,^o+/^ G(k,aj)e-''^^'. Using (g}, one has CT 



For simplicity, we consider first the case of uniform filling 
in the lattice, i.e. an integer number no of atoms per site, and 
we assume the system to be at zero temperature and in the 
insulating phase. In the limit of zero tunneling, the ground 
state wavefunction is a perfect MI, i.e. a product of number 
states at each site, and its Green function Go can be calcu- 
lated exactly (see appendix). The lowest-lying excited states 
of the system are "particle" and "hole" states, where a supple- 
mentary particle is added (respectively removed) at one lat- 
tice site. Creating these excitations costs a finite interaction 
energy, respectively = ;7no and = UUiq - 1) ill. 

To calculate the quasi-momentum distribution for a finite 
tunneling f, many-body techniques can be applied to obtain 
the single-particle Green function, G(k, oj). Using a path in- 
tegral approach, several authors 

O 111 El d El have 
been able to calculate the Green function of the Mott insulator 
within the RPA, 
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The poles E'j^^ of the Green function are the quasi-particle 
energies fl^ 
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- + f/(no- -) + -£)k[no], 
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In Eq. (|5Jl, fk - -2tYjv=x.v.z^'^^(^vd) is the dispersion re- 
lation for a free particle in the tight-binding limit, and 

Dk[no] - -yjtl^ + 4fkf/(no + 5) + U^- The particle weight is 

Zk - i^E^^ + ^) /^k[«o]- In the Appendix, we present an 
alternative derivation of ^ based on the equation of motion 
method, which follows closely Hubbard's method ilTl fT3l . 
Here we will simply comment on the physical picture behind 
this approach. The RPA considers that the particle/hole na- 
ture of the low-lying excitations is not significantly changed 
by introducing a finite tunneling (in technical terms, the self 
energy remains approximately the same as in the f ^ limit). 
The first effect of tunneling is to introduce a finite amount of 
particle/hole components in the the many-body ground state 
wavefunction. In the form given in 0], corresponding to a 
first order calculation, a particle/hole pair necessarily occu- 
pies two neighboring lattice sites due to the particular form 
of the tunneling hamiltonian. Through higher-order tunneling 
processes captured by the Green function (jS), the particle and 
the hole forming the pair can tunnel independently. As a re- 
sult, the pair acquires a mobility through the lattice, and may 
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'f + U{no + i) 



^tl + 4fkC/(no +k) + U^ ^ ^ 
To first order in tjU , this reduces to 

S{k) ^no- 2no(no + DtJU, 



(10) 
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also obtained in |8j|] by calculating the many-body wave func- 
tion perturb atively. We find that the two predictions rapidly 
converge. For example, they differ by less than 10 % for 
U/zt > 6.6, 11.6 and 16.9 for no -1,2 and 3, respectively. 
These values have to be compared to the respective critical 
values for MI formation, U/zt = 5.83,9.89, 13.93. This in- 
dicates that the coherence beyond nearest neighbors is rather 
rapidly lost as one goes further into the MI phase. However, 
the visibility itself remains finite in a substantial range of U/zt, 
implying a persistent short-range coherence. 



IV. COMPARISON WITH THE EXPERIMENTS 

To compare with the experiments reported in several 
features have to be taken into account. First, only the col- 
umn density n is accessible experimentally, i.e. the density 
integrated along the probe line-of-sight (which we take here 
parallel to the z axis). Second, the visibility is experimentally 
deduced from two points according to 



exp 
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To eliminate the Wannier envelope, nmax and nmin are mea- 
sured from two points at the same distance from the cloud 
center (see Fig. |2j, so that the envelope automatically can- 
cels out. For example, nmax is found at point {2n/d,Qi) and 
n^iin at {2n/ y/2d, 2n/ V2t/). This reduces the visibility com- 
pared to the usual definition. In the theoretical calculation it 
is straightforward to account for these two effects. 

The third effect, the shell structure of the MI, is handled 
here in an approximate way. In the numerical calculations, 
the shell distribution always includes small regions with non- 
integer filling, which the theory above cannot handle. How- 
ever, these domains are small, and have a strongly depleted 
superfluid component, so that we do not expect them to have 
a large effect on the visibility. Therefore, we approximate the 
density distribution by a "ziggurat"-like profile, where only 
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FIG. 2: Measurement of visibility. The interference pattern shown in 
the left graph corresponds to a lattice depth of 8 iiR, in the superfluid 
regime. The right graph indicates the geometry of the reciprocal lat- 
tice. Gray areas are the first and second Brillouin zones (projected 
in the image plane), and the white dot indicate the position of the 
maxima of the interference pattern. Along the circle, the Wannier 
function envelope takes the same value, and we measure the interfer- 
ence "minimum" at the intersection of this circle and of the diagonal 
of the lattice square lattice, indicated by the black dot. 



MI shells are present. The actual extension of each shell is 
calculated as if t were zero, taking the external potential into 
account I25I1 . In Fig. ^ we compare the profile in this ap- 
proximation (dotted line) with the numerically calculated one 
(solid line). For large lattice depths, both agree reasonably. 
Note that the density profile still depends weakly on the lat- 
tice depth through the external confinement [see Eq. 0]. 

The momentum distribution deduced from Eqs. jSllOt is 
averaged over the distribution of atoms to compare with the 
experimental data (see for details on the experiment). The 
results are plotted versus lattice depth in Fig.|3| for two differ- 
ent atom numbers in the lattice. For the lowest atom number 

= 2.2 X 10^, we calculate that only /iq = 1 and no - 2 shells 
are present. For the largest = 5.6 x 10^, a core with «o = 3 
atoms per site is also present. Note that in the latter case, the 
actual density distribution might deviate more from the cal- 
culated one, due to three-body losses in the «o - 3 region. 
We find that the calculation agrees with the measured visibil- 
ity within 20 % for Vq < 22 Er. The theory curves terminate 
when the MI shell with highest filling disappears, as it is re- 
placed by a large superfluid core not described by our theory. 
Note that the calculation does not include any free parameter. 

However, we consistently find that the calculated value lies 
below the measured visibility for large lattice depths Vq < 
22 Er. Moreover, the deviation increases with increasing lat- 
tice depths, which shows that the superfluid shells play little 
role in determining the visibility for such large lattice depths, 
as assumed in our calculation. In Fig. |3 the fractional de- 
viation of the calculated visibility from the measured one is 
plotted versus lattice depth for four data sets. Remarkably, al- 
though the atom numbers are rather different from one data set 
to another we find a common trend in the data. On the other 
hand, this observation also suggests that a breakdown of adia- 
baticity occurs for the particular ramp used in the experiments 
to increase the laser intensity to its final value, a point already 
identified in jsj]. We conclude that, perhaps surprisingly, the 
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FIG. 3: Comparison between the measured and the calculated vis- 
ibility. The upper and lower graphs correspond to atom numbers 
N = 2.2x and A' = 5.6x 10^ respectively. The dotted and dashed 
lines indicate the calculated visibility for homogeneous MI with fill- 
ing factor no = 1,2. The solid lines are calculations including the 
inhomogeneous shell distribution. Typical standard deviations for 
the experimental data are 1 % or below. Our calculation of the equi- 
librium distribution at zero temperature indicate that in case a, only 
MI regions with no = 1 and wo = 2 atoms per site form, whereas in 
case b, a core with ;7o = 3 is also present. 



visibility in the MI may be a sensitive probe of the many-body 
dynamics of the superfluid-to-insulator transition. 



V. CONCLUSION 

In conclusion, we have derived in this paper a theoretical 
expression for the interference pattern of a Mott insulator af- 
ter release from the optical lattice and a time of flight. Our 
calculations take deviations from perfect filling due to a finite 
tunneling into account, and use a simplified but realistic model 
of the shell structure of the MI. Good agreement with our ex- 
perimental data reported in 18] is found, at least for moderate 
lattice depths. For very large lattice depths, an increasing de- 
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FIG. 4: Fractional deviation of the calculated visibility from the mea- 
surements. Symbols denote different data sets with different atom 



numbers and densities (□: N = 2.2 x 10^, 
A' = 4.3x 10^•: A' = 5.9x 10^). 
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viation points to non-adiabatic effects in the conversion from 
a condensate to an insulating state, which could in principle 
be studied by the method presented here. Nevertheless, in 
view that no free parameter is included in the theory, we con- 
clude that the momentum distribution jlOt describes the sys- 
tem well. This supports the physical picture of the system as 
a (dilute) gas of partice/hole pairs, mobile through the lattice, 
on top of a regularly filled Mott insulator. Furthermore, the 
validity of the RPA to describe their behavior is qualitatively 
verified. 



Our calculation neglects entirely the superfluid component, 
which is correct only for large lattice depths where the sys- 
tem is almost completely insulating. Recently, several authors 
frMl26ll27ll have proposed to modify the standard mean-field 
description I19ll2lll to better account for long- and short-range 
coherence. It would be interesting to compare the predictions 
of those approaches with our data for lower lattice depths, 
where the system is expected to be a strongly depleted su- 
perfluid, and therefore amenable neither to a Bogoliubov-like 
description nor to a strongly interacting one as provided in 
this paper. Also, an investigation of finite temperature effects 
fllll would be useful. In particular, an interesting question is 
to know whether the visibility measurements presented here 
could be used for thermometry in the lattice. 



We would like to thank Dries van Oosten, Paolo Pedri 
and Luis Santos for useful discussions. Our work is sup- 
ported by the Deutsche Forschungsgemeinschaft (SPP1116), 
the European Union under a Marie-Curie Excellence Grant 
and AFOSR. EG acknowledges support from a Marie-Curie 
Fellowship of the European Union. 



APPENDIX A: GREEN FUNCTION OF THE 
HOMOGENEOUS MOTT INSULATOR IN THE RANDOM 
PHASE APPROXIMATION 

In this Appendix, we present a derivation of Eq. (|8} using 
the equation of motion approach. The single-particle Green 
function is defined at zero temperature as 



Gijit) = -i{Taim]it)} 



(Al) 



= -/0(f)<a,(O)a](/)> 
-ie{-t){a](t)aim, 

where T is the time-ordering operator and 9 is the Heaviside 
step function. Since we consider a time-independent and ho- 
mogeneous system, we take a Fourier transformation of this 
equation with respect to space and time (denoted by the sym- 
bol 'F), and define 



G{k,cj) = r[Gij{t)\ 



(A2) 



In the frequency-momentumrepresentation, the Heisenberg 
equation of motion ifij^Gij(t) - ^H, G,j(f)] takes the form 



(ha> + ju)G(k, oj) 



1 + fkG(k, w) 

-ir[{TnMO)a](t)}]. (A3) 



The last term on the right hand side of jA3> is usually rewritten 
as E(k, w)G(k, co), where E is the self-energy. This gives the 
expression 



G(k,w) = 



1 



Tico + /i - f k ~ 2(k, cj) 



(A4) 



Let us first assume that no tunneling is present (t - 0). In 
this case, the Green function can be obtained exactly from 
its definition JAll i and the ground state wave function \'^} = 
Yli \no}i, where each site is in the Fock state |/io). The result 
Go(ai) is independent of momentum, and reads 



Go(w) 



no + 1 



no 



Hcij + jJ. — Uno tioj + fi - U(no - 1) 



(A5) 



In this self-interacting Umit, we can rewrite Eq. JA5> as 
Go(w)"' = + 1-1- So(w), with the self energy 



So(w) = 2Uno 



U\)(no + 1) 
hco + jj. + U 



(A6) 



This expression is exact in the f — > limit, and coincides with 
the one found in fl^ . The first term is simply the Hartree- 
Fock energy per particle for uncondensed atoms (hence the 
factor of 2), whereas the second term - which has the same 
order of magnitude at low energy - accounts for the coiTela- 
tions between particle that drive the system into the perfectly 
ordered ground state. 

If we now restore a finite tunneling, but still consider a sys- 
tem in the insulating phase, a reasonable approximation is to 
assume that the self energy is not changed with respect to the 
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strongly interacting limit. We comment on this approximation which has a typical RPA form. Using the explicit result for 
in the test. Making this approximation yields Go, we obtain after some algebra Eq. (jSJl in the text, which 



1 

G(k, w) 



1 -f(k)GoM' 



(A7) 



explicitly displays particle and hole components. 
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